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ABSTRACT 

We develop a new method to identify YSOs from star-forming regions using the photometry data 
from Spitzer's c2d Legacy Project. The aim is to obtain YSO lists as complete as possible for study- 
ing the statistical properties, such as Star Formation Rate (SFR) and lifetimes of YSOs in different 
evolutionary stages. The largest obstacle for identifying YSOs comes from background galaxies with 
similar SEDs to YSOs. Traditionally, selected color-color and color-magnitude criteria are used to 
separate YSOs and galaxies. However, since there is no obvious boundary between YSOs and galaxies 
in Color-Color Diagrams (CCDs) and Color-Magnitude Diagrams (CMDs), those criteria may exclude 
faint YSOs near the boundary. In this paper, we separate the YSOs and galaxies in multi-dimensional 
(Multi-D) magnitude space, which is equivalent to using all variations of CMDs simultaneously. Com- 
paring sources from molecular clouds to Spitzer's SWIRE data, which have negligible amount of 
YSOs, we can naturally identify YSO candidates locating outside of the galaxy populated regions 
in the Multi-D space. In the five c 2d-surveyed c l ouds, we select 322 new YSO candidates (YSOc), 
miss/exclude 33 YSOc compared to Evans et al. (2009) and result in 1313 YSOc in total. As a re- 
sult, SFR increases 28% correspondingly, bu t the lifetimes of YSOs in different evolutionary stages 
remain unchanged. Comparing to theories (Krumholz & McKee 2005), our derived SFR suggests 
that star formation in large scale is dominated by supersonic turbulence rather than magnetic fields. 
Furthermore, we identify 7 new Very Low Luminosity Objects. 

Subject headings: ISM: clouds - Stars: luminosity function, mass function - Stars: pre-main sequence 
- Stars: protostars - brown dwarfs - Stars: formation 



1. INTRODUCTION 

A full census of the Young Stellar Object (YSO) pop- 
ulation in star-forming regions is essential for accurately 
calculating statistical properties, such as Star Formation 
Rate (SFR) and lifetimes of YSOs in different evolution- 
ary stages, which are fundamental parameters for assess- 
ing the physical mechanism of global star formation. The 
idea that SFR shoul d be related to gas surface density 
was first proposed by Schmidt ( 1958 ) and the rel ation is 
measured in a galaxy sample by ]Kennicutt ( 1998 ) , which 
is known as Kennicut -Schmidt relation. To examine this 
relation in our Galaxy, Evans et al.| ( |2009"| count the num- 
bers of YSOs in Giant Molecular CTouds (GMCs), and 
assume a mean mass and formation timescale for these 
YSOs. They find that their measured SFRs are higher 
than that indicated by the Kennicut-Schmidt relation, 
and claim that is because the Kennicut-Schmidt rela- 
tion applies to av erage over much l a rger r egions than 



not been able to determine which is the domi nated sup 
porting force. Krumholz & McKee ( 2005 ) and Krumholz 



individual clouds. Heiderman et al. (2010) extend the 



sample and obtain a simil ar resultto t hat fr om Evans 



etaL|j2009j. In addition, |Lada et al.j ( |2010[ ) find that 
the SFR is linearly proportional to the mass of cloud 
above an extinction threshold and this relation shows 
excellent agreement between galactic and extragalactic 
star-forming activity. 

However, what physical mechanism determines the 
SFR is still unclear. It has long been known that the 
observed SFR is too low for a GMC collapsing on its 
free-fall timescale with out support against gravity ( Zuck- 
erman & Evans 1974). Supporting force such as super- 
somc turbulence and/or magnetic fields has been used to 



& Tan (2007) derive analytical relations between SIR 
and viral parameter (the ratio of kinetic energy and grav- 
itational energy), and show that SFR is linearly propor- 
tional to viral parameter for magnetic field dominated 
scenario while SFR decreases exponentially versus viral 
parameter for turbulence dominated scenario. Here we 
attempt to obtain SFR as accurate as possible by devel- 
oping a new YSO census method. In addition, a better 
census method will provide better estimates in the life- 
times of YSOs in different evolutionary stages, which are 
usually estimated by comparing the fractions of YSOs 
in each stage, and successful theories should be able to 

explain the observed YSO lifetime scales. 

Spitzer Space Telescope ( Werner et al.||2004 ) provides 
high sensitivity surveys with the Infrared Array Cam- 
era ( IRAQ at 3.6, 4.5, 5.8, and 8.0 ^m ( |Fazio et aT 
20041 and the Multib and Imaging Pho tometer (MIPS) 
70, and 160 /jm ( |Rieke et al.|2004| , which allow us 
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to search for YSOs in star-forming regions. Identifying 
YSOs is usually achieved by removing stars and back 



ground galaxies in the data (Young et al. 2005 



et al.||2007a[ iRebull et all 2007 
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Stars can 
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be eas- 



ily selected by fitting SEDs with reddened stellar atmo- 
sphere. However, the SED morphologies of galaxies are 
very similar to those of YSOs at infrared wavelengths 



(Harvey et al. 



2006J, and Spitzer can detect substan- 
tial amount of' background galaxies because of its high 
sensitivity. Therefore, separating background galaxies 
and YSOs becomes a difficult problem. Several methods 
explain thelow" SFR^howeTer," otee^t^aTd^tl hlw have been develo P ed to identif y YSOs from the molecular 



clouds (Harvey et al.||2007a Gutermuth et al 



2005 Re 



bull et al.|2010[) using Spitzer data. These methods iden- 
tify YSOs by eliminating the background galaxies based 
on comparing the distributions of the observed data with 
galaxy data in Color-Color Diagrams (CCDs) and Color- 
Magnitude Diagrams (CMDs). However, there are no 
obvious boundaries between YSOs and galaxies in CCDs 
and CMDs. Hence, based on different considerations, dif- 
ferent works use different sets of CCDs and CMDs, and 
set their own boundaries in CCDs and CMDs to identify 
YSOs. Therefore, sources located close to the boundaries 
are possibly classified as different objects using different 
methods. 

A large YSO survey toward five nearby molecular 
clouds has been done by the Spitzer Legacy Project 
"From Molecular Cores to Planet Forming Disk" (c2d; 
Evans et al. 2003). Th e c2d project used a n unnormal- 
ized galaxy probability ( |Harvey et al.|200"7a I to eliminate 
possible galaxies and identity YSO candidates (YSOc) 
from the five clouds. (In the paper, YSOc is used for 
sources selected by our method or the c2d project, which 
may or may not be confirmed as true YSOs by other 
studies.) This galaxy probability is primarily calculated 
from the locations of sources in three CMDs by com- 
paring to the Spitzer Wide-Area Infrared E xtragalactic 
Survey (SWIRE) data ( |Lonsdale et al.|2003[ ), which con- 
tain negligible amount of YSOs because the observations 
were towards high galactic latitude. Since this galaxy 
probability was calculated only for sources detected in 
all IRAC and MIPSI (24 (mi), the c2d YSO catalog 
may miss faint YSOs which are undetected at one or 
more of the five bands. In this paper, we develop a new 
YSO identification method and apply it to the five c2d- 
surveyed molecular clouds, providing a more complete 
YSO catalog in these clouds. 

With a relatively complete YSO sample, we will be 
able to identify more Very Low Luminosity Objects (VeL- 
LOs), which are defined as protostars with internal lu- 
minosity, L; nt , smaller than 0.1 Lf?>. The first VeLLO, 



L1014-IRS was o bserved by |Young et al | ( |2004[ ). pun- 
ham et ah] ( |2008[ ) then design a set of color criteria and 
also use the c2d galaxy probability to select faint Class 
and early Class I sources as VeLLO candidates. The 
faint nature of VeLLOs suggests that they could be ei- 
ther very young protostars or very low mass protostars, 
but recent works have shown that the low luminosity 
can be explained by protosta rs in the quiescent pha se of 
episodic accretion proce sses ([Dunham et al1|2010[ ). To 
identify more VeLLOs, Dunham et al.| ( |2U08p use the 
galaxy probability from c2d project but adopt a higher 
cutoff, which allows to find more faint sources; the galaxy 
probability are calculat ed based on the general fact that 



th e galax ies are faint (Dunham et al. 2008 Harvey et 
al. 2007a). However, the change ot cutoff may result in 



higher galaxy contamination and more analysis are es- 
sential for confirming those VeLLOs of |Dunham et al.| 
(2008). Therefore, a method that identifies faint YSOc 
naturally may help to reduce the galaxy contamination. 

We develop a new method, multi-dimensional (Multi- 
D) method to identify YSOs naturally. While a CMD 
(magA, magB-magC) contains only information from the 
three consisting magnitudes, Multi-D magnitude space 
contains all the information of the consisting magnitudes. 



We treat the number distribution of our galaxy sample 
in Multi-D magnitude space as galaxy probability distri- 
bution and sources located in the region without galaxies 
are classified as YSOc. We apply this method to the five 
c2d-surveyed clouds and the number of YSOc sample is 
increased by a factor of 28%. Therefore, the SFR we cal- 
cul ated in this work ar e 28% larger than that calculated 
by Evans et al.| ( 2009 ) due to the increased number of 
YSOc, but the life times of YSOs in different evolution- 
ary stages remain unchanged due to the similar ratios 
between different stages of newly found YSOc. We de- 
scribe the data we used in §2 and our YSO identification 
method in §3. In §4, we discuss how reliable our new 
YSOc are. The statistical analysis and discussion are 
presented in §5 and the conclusions are summarized in 
§6. 

2. DATA 

The data we used in this paper are from two Spitzer 
Legacy Projects, c2d and SWIRE. c2d project observed 
five nearby molecular clouds and SWIRE project ob- 
served extragalactic fields that contain negligible amount 
of YSOs. Our aim is to identify YSOs from the c2d data 
by comparing them to SWIRE data in Multi-D magni- 
tude space. 

2.1. c2d Data 

Perseus, Serpens, Ophiuchus, Lupus (I, III and IV) 
and Chamaeleon II have been observed by c2d project in 
IRAC1-4 (3.6, 4.5, 5.8 and 8.0 fim) and MIPS1-3 (24, 70 
and 160 /im) bands. Hereafter, we use IR and MP to rep- 
resent IRAC and MIPS, e.g., IR1 as IRAC1. All data are 
processed with the c2d standard pipeline and these pro- 
cesses ar e described in the c2d data delivery document 
in detail ( |Evans et al.|2007 ). The pipeline extracts point 
sources from images at all bands and performs photo- 
metric measurements. The Spitzer photometry data are 
merged with 2 Micron All-Sky Survey (2MASS) data (J, 
H, K bands) into source catalogs using a 2" matching 
radius. These source catalogs and images are released 
to Spitzer Heritage Archive. Detail studies of individual 
clouds using these catalogs have been presented in several 
papers (e.g., Perseus: J0rgensen et al. 2006; Serpens: 
Harvey et al. 2007a; Ophiuchus: Padgett et al. 2008; 
Lupus: Merin et al. 2008; Chamaeleon II: Porras et al. 
2007). Two kinds of catalogs, FULL and HREL (High 
Reliability), are provided by the c2d project. FULL cat- 
alogs include all detected sources but some of them could 
be fake detection misidentified by automated photometry 
processes. Sources in HREL catalogs have to be detected 
in one band with S/N > 7 accompanied by a detection 
in another band with S/N > 5. In this paper, we select 
YSOc only from HREL catalogs in order to avoid fake 
sources. 

We adopt the photometry results from the c2d cat- 
alogs for our analysis. Source flux in the c2d catalogs 
are measured with Point Spread Function (PSF) fitting. 
An important factor, image type (imtype), is used to 
indicate the quality of PSF fitting. Sources can be well 
fitted with PSF are assigned imtype = 1 and sources 
with imtype =^ 1 have high uncertainty in photometry 
resulting from low S/N ratio, extende d size, or non- 
circular shape, etc (lEvans et all 120071). For example, 
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Fig. 1. — The YSO identification process. Sources is consistent 
with the criteria in the block will enter the block connected by "Y" 
arrow and "N" for inconsistent. 

imtypc = 7 indicates that the S/N ratio is too low 
to be tested for shape and we found that it is mostly 
due to the confusion from surrounding cloud structure. 
Therefore, we will check the images of sources with 
imtype ^1 to determine whether they are real sources 
(see §3). 



2.2. SWIRE Data 
The SWIRE data used in this work are the same data 



used by Harvey et al. (2007a), which was processed with 
c2d pipeline to eliminate the systematic uncertainty in- 
troduced from different data reduction processes. There- 
fore, it can be used to compare with the c2d data fairly. 
The SWIRE data used here are the observation toward 
ELAIS Nl cxtragalactic field which is at high latitudes 
(b ~ 44°) and the observed field has coverage of 5.3 deg 2 
using both IRAC and MIPS. Since the YSOs are concen- 
trated in the galactic plane, there should be almost no 
YSOs in the SWIRE region. In order to acquire a galaxy 
sample as pure as possible, we removed the stars in the 
SWIRE data using the same method to that for c2d data 
(§3.2). Thus, it can be used to represent a full collection 
of the SED features of background galaxies. 

3. YSO IDENTIFICATION 

3.1. Selection Process 

Our YSO identification procedure is shown in Fig. 1 
and described below. Basically, this process identifies 
YSOc by eliminating the non-YSO sources. 



1. The main sequence stars were first removed by SED 
fitting. The fitting process is the same as that in |Har-| 
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. Oliveira et al. (2009) 
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Fig. 2. — The giant star selection criteria. Green points are 
the SWIRE data exc luding stars and red points are giant stars in 
Serpens identified by |Qliveira et a l. (2009). Sources located in the 
region below the dashed line are classified as giant stars. 



vey et al.| ( |2007b[ ) and |Evans et al.| ( |2007[ ), which se- 
lect stars with reddened stellar atmosphere SE D tem- 
plates. The extinc t ion law with Ry = 5.5 from|Wcin- 
gartner & Draine (2001) is used to obtain the best 
fit extinction value, Ay- Weingartner and Draine Ry 
= 5.5 model is suggested to be a good description 
for d usts in dense molecular clouds (Chapman et al. 



2009) 



Since our galaxy sample (SWIRE) is observed toward 
regions with negligible extinctions, we de-reddened the 
whole c2d catalog in order to compare it with our 
galaxy sample. We use the sa me computer progr am 



developed by the c2d project (Evans et al. 20071 to 
construct the extinction maps from trie extinction of 
background stars. The average errors of all pixels in 
each extinction map are 0.03-0.05 magnitude for each 
cloud. Sources in the c2d catalog are then de-reddened 
according to the extinction value of its position in the 
map. We use the de-reddened fluxes in this paper ex- 
cept for the analysis with a (section 5.3.2). 



3. Giant stars were removed using a CCD, IR2 
versus IR3 - MPl (Fig. 2), as described in §3.3. 



IR3 



We found that some sources labeled as "U" (unde- 
tected) for MPl in the c2d catalog are in fact satu- 
rated. Since bright MPl flux will make the source far 
away from the galaxy populated region in the Multi-D 
array, we believe these sources are very young YSOs 
rather than galaxies. Therefore, we classify these 
MPl saturated sources as YSOc. To find the satu- 
rated sources, we search all sources labeled as "U" at 
MPl with pixel value larger than 800 MJy sr _1 within 
7.66"radius. The cutoff 800 MJy sr _1 is determined 
from experience by checking the pixel values around 
many saturated sources, and we found that it is small 
enough for us to pick up all saturated sources and also 
not too small for us to obtain too many fake sources. 
The 7.66" radius is the 3tr radius of MPl PSF. Then, 
we examine the images of selected sources by eyes to 
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confirm whether they are indeed saturated sources and 
the excess bright flux is indeed from a saturated source 
but not from a bad pixel. The typical images of sat- 
urated sources are PSF-like with saturated holes (no 
values or low values) in centers and bright rings around 
the centers. 

5. To eliminate background galaxies, many Multi-D ar- 
rays are constructed for the SWIRE data to calculate 
the smoothed galaxy density. According to the loca- 
tion of the c2d sources in the Multi-D arrays, we cal- 
culate four galaxy probabilities: Gal prob 1, Gal prob 
2, Gal prob lp and Gal prob 2p, where Gal prob 1 and 
Gal prob 2 are the galaxy density in (J, Ks, IR2, IR4 
and MP1) and (IR1, IR2, IR3, IR4 and MP1) arrays 
or their subarrays (see §3.2 for detail), and Gal prob 
lp and Gal prob 2p are the galaxy density from the ar- 
rays discounting the bands with imtype ^ 1. Sources 
with both Gal prob 1 and 2 > 1 are classified as galaxy 
candidates and removed. 

6. The remaining sources are with Gal prob 1 < 1 or 
Gal prob 2 < 1, thus are potential YSOc. Images of 
sources with Gal prob 1 (2) < 1 but Gal prob lp (2p) 
> 1 (meaning an unreliable detection would affect its 
galay probability) will be examined by eyes. Sources 
with unreliable photometry such as jet knots or cloud 
structures are removed (see §3.2). 

7. For all the remaining sources, we also check the IR1 
images to eliminate the nearby resolved galaxies which 
can be bright and thus with low galaxy probability. 
Eleven nearby galaxies are removed through this pro- 
cess. 

8. The remaining sources are classified as YSOc. 

3.2. Multi-D Array Construction and Galaxy 
Probability Calculation 

The heart and soul of our new method is to use the 
galaxy density in the Multi-D space as an indicator of 
the "galaxy probability" ; anything outside of the galaxy 
populated region would have galaxy density smaller than 
1 and thus will be selected as YSOc. We assume that 
SWIRE data is large enough to contain a complete set of 
galaxy sample. Note that our "galaxy probability" here 
is not a real probability for a source to be a galaxy, it 
is in fact an unnormalized number which indicates how 
many galaxies near a specific position in Multi-D space. 
In this section, we describe the details for constructing 
Multi-D space and defining the galaxy probabilities. 

The data we used contain ten bands (J, H, Ks, IR1, 
IR2, IR3, IR4, MP1, MP2 and MP3). Although a high- 
dimensional array includes more information than that 
of a low-dimensional array, the maximum dimension of 
the data array we can handle is limited by our computing 
resource. If we want to calculate the galaxy probability 
from all 10-band data using a 10-D array, our computer 
will need random access memory of ~10 8 G byte (with 
the cell sizes we choose, see below), which largely exceed 
the limitation of our computer. Therefore, we choose to 
analyze our data with two 5-D main arrays construct- 
ing from the bands that can represent the main SED 
features; one contains J, Ks, IR2, IR4 and MP1 bands 
(corresponding to Gal prob 1) and the other contains 



all IR1-4 and MP1 bands (corresponding to Gal prob 2). 
MP2 and MP3 bands are excluded due to their low detec- 
tion numbers and their poor angular resolutions (18" and 
40", respectively). The J-MP1 array covers the widest 
range of wavelengths and the IR1-MP1 array is designed 
for selecting YSOc that have no 2MASS detections. The 
IR1-MP1 array is necessary for identifying young embed- 
ded YSOc since they are very faint at short wavelengths 
and are often undetected in 2MASS bands. 

Because the galaxy locations are discrete in the Multi- 
D arrays, we need to carefully choose the cell size and 
also employ a smoothing process to produce a smoothed 
galaxy distribution, so that the galaxies fall in between 
SWIRE sources will not be classified as YSOc. We find 
that using cell size of 0.2 mag for all bands is adequate for 
our data. We smooth the data using Multi-D Gaussian 
beams with a peak value of 1, a standard deviation a = 
2 cells (0.4 magnitude) in every dimension and an outer 
cutoff at radius of 7 cells. The galaxy probability for a 
source located at (xi, X2, X3, X4, X5) is defined as the total 
number of galaxies, A(xi, X2, X3, X4, X5), after smoothed 
with Multi-D Gaussian beams, 

xi=x t +7 x 2 =x 2 -7 x 3 =x 3 +7 x 4 = x 4 + 7 x^=x 5 + 7 

E E E E E 

x[=x 1 -7 x 2 =x 2 -7 13 = 13-! x^=x 4 -7 x^=x s -7 



(xi - x[) 2 + (x 2 - x' 2 ) 2 + (x 3 - x' 3 ) 2 + (xi - x 4 ) 2 + (x 5 - x' s ) 2 < 49. 

(1) 

Therefore, Gal prob 1 = A smooi /, (J, Ks, IR2, IR4 and 
MP1) and Gal prob 2 = A smooth (IR1, IR2, IR3, IR4 
and MP1). If a galaxy is surrounded by other galaxies, 
the galaxy' galaxy probability should be larger than the 
peak value of smoothing beam, which is 1, because the 
smoothing process will accumulate the number of galax- 
ies from nearby cells. Therefore, the surface with galaxy 
probability — 1 is expected to enclose the galaxy popu- 
lated region in Multi-D magnitude space except isolated 
galaxies. For an isolated galaxy, its galaxy probability 
will be one. Out of 135400 SWIRE galaxies, there are 
only three galaxies with Gal prob 1 equal to one and no 
galaxies with Gal prob 2 equal to one, which indicates 
that the galaxy sample is well smoothed in the Multi-D 
arrays. Therefore, sources located out of the smoothed 
galaxy distribution will have galaxy probability < 1 and 
can be classified as YSOc. 

Since a large number of sources are not detected in 
all five bands of these two 5-D arrays, for these sources, 
Gal prob 1 and Gal prob 2 are calculated from 3- and 
4-D subarrays instead. Thus, sources with detections in 
three to five bands of the five bands for Gal prob 1 arc 
assigned Gal prob 1 calculated from 3- to 5-D magnitude 
arrays, respectively, and same as Gal prob 2 (the detec- 
tion threshold is set to S/N=2). As a result, 15 subarrays 
(five 4-D and ten 3-D subarrays) for each 5-D array are 
constructed. The array cell size and/or the smoothing 
lengths of the subarrays need to be modified, because 
the galaxy distribution is condensed from 5-D to 3- or 
4-D and thus the number of galaxies in 3-D or 4-D array 
cell are much larger than that in 5-D array cell. How- 
ever, it is very difficult to make the galaxy probability 
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Fig. 3. — Reducing beam size from 2-D to 1-D. The two circles 
indicate the beam size in 2-D and the two half-ellipses indicate the 
projected beam size in 1-D. The "galaxy probability" along the 
dashed line, diagonal direction, in 2-D will be the same to that 
along x-axis in 1-D if a reducing beam length of ^/l/2<7 is applied. 

of each source be the same in different dimensions and 
fairly compare them in different dimensions. Therefore, 
we modified the smoothing beam to make the threshold 
of galaxy probability, 1, a reasonable number for sep- 
arating YSOc and galaxies in arrays with different di- 
mensions. Here we demonstrate the principle of beam 
size modification in the case of 2-D to f-D. We assume 
that there are very few galaxies close to the surface of 
galaxy populated region where the galaxy probability is 
~1. Considering two galaxies locating near the surface 
and along the diagonal direction, we find that the galaxy 
probability profile in 1-D is exactly the same as that in 
2-D along the diagonal direction if we reduce the stan- 
dard deviation of Gaussian beam in 2-D by a factor of 
(Fig. 3). 

Therefore, we use the standard deviation a of -\/3/5 x 
2 cells (the o~ in 5-D array is 2 cells) for the 3-D arrays and 

•\/4/5 x 2 cells for the 4-D arrays, and we keep the cell 
size unchanged. Such modification will produce similar 
galaxy probability distribution in the regions close to the 
boundary of galaxy populated regions for arrays with 
different dimensions, thus the same threshold of galaxy 
probability (i.e., 1) can be used to identify YSOc. 

We found that a lot of sources with PSF fitting flag 
"imtype ^ 1" (§2.1) have contaminated flux measure- 
ments due to jet knots, nearby bright sources and ex- 
tended cloud structure, which may change the locations 
of the sources in the Multi-D arrays. Thus, we define Gal 
prob lp and Gal prob 2p which are calculated from sub- 
arrays only using the detections with imtype = 1 (well- 
fitted as a point source), and use them along with Gal 
prob 1 and Gal prob 2 for separating YSOs and galax- 
ies. Since imtype=l indicates that the source has a more 
precise flux measurement, source with Gal prob lp (or 
2p) < 1 implies that this source locates out of the galaxy 
populated region in Multi-D space with high confidence. 
If a source has Gal prob 1 (or 2) < 1 but Gal prob lp 
(or 2p) > 1, it is obvious that the bands with imtype ^ 



1 make the sources deviated from the galaxy populated 
region. For these sources, images of IR1-IR4 and MP1 
bands are examined with eyes. If a reliable point source 
is found in the images at bands with imtype ^ 1, the 
Gal prob 1 (or 2) calculated containing such bands is re- 
liable and the sources are classified as a YSOc. A large 
amount of examined sources were found to have unre- 
liable photometry measurement (70%-90% for different 
clouds) due to jet knots, nearby bright sources and ex- 
tended cloud structure, and are removed from our YSOc 
list. 

3.3. Giant Star Contamination 

Giant stars could be bright and have similar SE Ds to 
Class II/III YSO in infrared ( |01iveira et aLp 009), thus 
are very difficult to be removed only using 1KAC and 
MIPS data. We are not able to exclude the giant stars 
by finding the natural boundary of them in the Multi-D 
array, because the giant star sample is too small to repre- 
sent a full collection of the SED features. From SWIRE 
data, we found that there are 30 bright sources that are 
not group with the majority of the SWIRE sources in 
the (IR2 - IR3, IR3 - MP1) CCD. Because we think that 
SWIRE data only contain galaxies, main-sequence stars, 
and giant stars, we suspect that these 30 bright sources 
are giant stars. Therefore, we define 3 criteria to exclude 
the giant stars from the c2d data: IR3 - MP1 < 2, IR2 
- IR3 > IR3 - MP1 and IR2 - IR3 > (Fig. 2). The 
numbers of giant stars selected from these criteria are 9, 
8, 5, 2 and 1 in Serpens, Ophiuchus, Lupus, Perseus and 
Chamaeleon II, respectively. However, these criteria do 
not remove 17 sources in Serpens which have b een iden- 

~" <|2009[ > based on 
ions of these 17 



tificd to be giant stars by Oliveira et al. 
the optical spectrum. Because the ioca 
sources in Fig. 2 are closer to the SWIRE sources than 
our selected giant stars, we suspect that these spectro- 
scopically identified giant stars may have small disk or 
dusty envelopes that increase their IR3 - MP1 color. We 
later removed these 17 sources from our YSOc list. Since 
we did not find any spectroscopic data for giant stars in 
Perseus, Ophiuchus, Lupus and Chamaeleon II, we are 
not able to remove spectroscopic identified giant stars in 
these clouds. However, because Serpens molecular cloud 
is located in the direction closed to the Galactic plane 
while other clouds are not, Serpens has the most serious 
giant star contamination among these clouds. Thus, the 
giant star contamination should be much smaller in the 
other four clouds. 

3.4. Advantages of Multi-D Method 

Using the Multi-D space as a whole instead of using 
limited numbers of selected CMDs to identify YSOs has 
the following advantages: 

1. Avoidance of selecting specific CMDs. In order to 
identify YSOs in molecular clouds, previous works 
adopted different CMDs as their selection criteria. Al- 
though these CMDs are all selected with justifiable 
reasons, it is difficult to argue that which CMD sepa- 
rates YSOs and galaxies best and how many CMDs are 
sufficient for the identification. The various selections 
of CMDs in use make the YSO identification results 
different by different works. Since Multi-D magnitude 
space includes the information from all the possible 
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Fig. 4. — All variations of MMDs consisting any two bands of IR1, IR2, IR3, IR4, and MP1 for Perseus. Red, blue and green points 
indicate the "Same" YSOc, newly identified YSOc and SWIRE galaxies, respectively. The arrows represent the diagonal direction which is 
related to source brightness. 
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Fig. 5. — The same to Fig. 4 for Serpens. 
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Fig. 6. — The same to Fig. 4 for Ophiuchus. 
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Fig. 7. — The same to Fig. 4 for Lupus. 
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Fig. 8. — The same to Fig. 4 for Chamaeleon II. 



TABLE 1 

Numbers of selected YSOc and the comparison with the c2d lists 



Cloud 


YSOc selected with 
Multi-D Method 
(Group A) 


YSOc in the c2d 
lists (Group B) 


YSOc in both 
Group A and B 
("Same YSOc") a 


YSOc in Group A, 
but not in Group 

B a 


YSOc in Group B, 
but not in Group A 


Perseus 


469 


385 


369 (78%) 


100 (22%) 


16 


Serpens 


296 


227 


218 (74%) 


78 (26%) 


9 


Ophiuchus 


367 


392 


287 (78%) 


80 (22%) 


5 


Lupus 


143 


91 


91 (64%) 


52 (36%) 


3 


Chamaeleon II 


38 


26 


26 (68%) 


12 (32%) 





Total 


1313 


1024 


991 (75%) 


322 (25%) 


33 



Note. — The YSOc number includes the three well-known sources that are not selected from Multi-D Method (Table 3). 
a Percentage is calculated with respect to the YSOc number in Group A. 



variations of the CMDs, we have no need to select the 
specific CMDs. Therefore, Multi-D method is rela- 
tively complete and vmbias compared to CMD meth- 
ods. 

2. A natural boundary of galaxy populated region in 
Multi-D space. We use the natural boundary con- 
structed from number distribution of our galaxy sam- 
ple. This boundary encloses a relatively accurate 
galaxy SED ensemble, which reduces the bias in iden- 
tifying YSOc near the boundary. 

3. Uncovering YSOs that cannot be found with all pos- 
sible CMDs. A CMD is equal or equivalent to a pro- 
jection in Multi-D magnitude space. However, even 
all the possible variations of the CMDs are used, it is 
possible that some YSOs may not be found if they are 
located in a region without galaxies in Multi-D magni- 
tude space but immerses in galaxy populated regions 
in all CMDs; that is, some YSOs cannot be revealed 



in any projections from Multi-D magnitude space (see 
Fig. 4-9). 

Therefore, Multi-D method provides the best oppor- 
tunity to identify YSOs as complete as a photometric 
dataset can offer. 

4. RESULT: NEW YSOC LIST 

We present our YSOc lists for the five c2d-surveyed 
clouds selected with the Multi-D method. The total 
number of our YSOc and the difference between our 
YSOc numbers and the c2d numbers are listed in Tabic 

I, which shows that we increase the total YSOc num- 
ber by 28% (from 1024 to 1313). The YSOc we selected 
in Perseus, Serpens, Ophiuchus, Lupus, and Chamaeleon 

II, along with their galaxy probabilities are listed in Ta- 
ble 2. Note that the smaller a galaxy probability is, the 
more distant the sources is located away from the galaxy 
populated region. 
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TABLE 2 

YSO CANDIDATES 



R.A. 


Dec 


Gal prob l a 


Gal prob 2 a 


Gal prob lp a 


Gal prob 2p a 


Ste den 


YSO prob 


c2d classification 15 


degree 


degree 


log(num) 


log(num) 


log(num) 


log(num) 


log(pc- 2 ) 


log(num) 




Perseus 


51.1968071 


30.4656086 


-1.99 


2.44 


-1.99 


2.44 


-0.06 


1.03 


- 


51.3304656 


30.7577636 


_ 


-0.06 


— 


-0.36 


0.47 


0.46 


— 


51.3313375 


30.5733784 


-1.36 


-co 


— 


- 


0.43 


0.95 


YSOc_star+dust(IRl) 


51.3430184 


30.7538668 


-oc 


-■oc 


- 


- 


0.49 


0.00 


YSOc_red 


51.4009155 


30.7543631 


-co 


-co 






0.44 


0.23 


red 


51.4020445 


30.7561528 


-co 


-co 


-'OC 


-co 


0.43 


0.00 


YSOc_red 


51.4092048 


30.4527241 


-0.46 


1.14 


0.38 


1.86 


0.04 


0.01 




51.4117940 


30.7350532 


-co 


-co 






0.51 


0.03 


YSOcjed 


51.4130155 


30.7328223 


-co 


-co 




-co 


0.52 


0.01 


YSOc_red 


51.6561277 


30.2578016 


-co 


-co 






0.17 


0.06 


YSOc_red 


51.8292253 


30.2883997 


-0.32 


-0.44 






0.63 


0.14 




51.8393548 


30.8312253 


-2.22 


-■oc 






-0.16 


0.00 




51.9093741 


30.2329495 


-co 


-co 






0.87 


0.87 


YSOc_star+dust(IR4) 


51.9109468 


30.2267923 


-1.37 


-0.31 


-1.19 


0.47 


0.91 


1.50 




51.9117155 


30.2235869 


-1.24 


-0.24 


-0.34 


0.81 


0.88 


1.09 




51.9128437 


30.2175414 


-co 


-co 


-co 


-co 


0.89 


0.00 


YSOc_red 


51.9228069 


30.3379885 


-oc 


-'OC 






0.43 


1.64 


YSOc_star+dust(IR2) 


51.9301209 


30.2080268 


-co 


-co 


-co 


-co 


1.03 


0.00 


YSOcjred 


51.9486420 


30.2012577 


-co 


-'OC 


-co 




1.18 


0.62 


YSOc_star+dust(IRl) 


52.0003793 


30.1463963 


-oc 


-'OC 






0.87 


1.36 


YSOc_star+dust(IR2) 



a We calculate four galaxy probabilities: Gal prob 1, Gal prob 2, Gal prob lp and Gal prob 2p, where Gal prob 1 and Gal prob 2 are the galaxy 
density in (J, Ks, IR2, IR4 and MP1) and (IR1, IR2, IR3, IR4 and MP1) arrays or their subarrays (see §3.2 for detail), and Gal prob lp and Gal 
prob 2p arc the galaxy density from the arrays discounting the bands with imtypc ^ 1. 

If the c2d source type does not contains "YSO" , the source is added to the YSOc list by [Evans et ahl ]2009[ based on ancillary data from 
literature. Sources without c2d classification, "— " , are newly identified YSOc from this work. 



We discuss the properties of our newly selected YSOc 
in §4.1 and argue that they are likely to be real YSOs. 
In §4.2, we present the reasons for excluding some YSOc 
from the c2d list and retrieve three missed established 
YSOs. In §4.3, we further discuss the YSOc with "PAH- 
emission" feature and label them as less reliable YSOc, 
which composes less than 3% of all YSOc. Finally in 
§4.4, we discuss the uncertainty resulting from our band 
selections for constructing two 5-D arrays. 

4.1. Analysis of The Newly Identified YSOc 

In order to examine whether our newly identified YSOc 
are likely to be real YSOs, we compare our list and the 
c2d list in the following three aspects. First, we show the 
distributions of YSOc from the two lists are consistent in 
all MMDs and the two lists of YSOc have similar lumi- 
nosity function. Second, we show that our YSOc stellar 
surface density are more consistent with the c2d's YSOc 
stellar density than galaxy surface density. Finally, we 
show that the distributions of the two YSOc lists are also 
consistent in the multi-dimentional color space. Note 
that these analyses can only give us the hints of that 
whether the newly identified YSOc are likely to be real 
YSOs or not. Without other data such as spectroscopy, 
we are not able to confirm whether they are real YSOs. 

4.1.1. Newly Identified YSOc in The Magnitude-Magnitude 
Diagrams and Their Infrared Luminosity Function 

Because we are not able to plot a Multi-D figure with 
more than 3 dimensions, we display the distribution of 
our newly identified YSOc in all possible MMDs, in order 
to examine whether their distributions in the MMDs are 
similar to the previously identified YSOc. For the same 
type of sources with similar SEDs, their locations in the 
CCDs will be close to each other and their locations in 
MMDs will be along the diagonal direction since the only 



difference between sources is the brightness (distance). 
Therefore, if the newly identified YSOc are real YSOs, 
they should align with the previously identified YSOc in 
the diagonal direction. Fig. 4-8 show the populations 
of newly identified YSOc (blue), galaxy sample (green) 
and YSOc identified by both c2d and our method (red) 
(hereafter we call these YSOc as "Same" YSOc) in all 
variations of MMDs consisting any two bands of IR1, 
IR2, IR3, IR4, and MP1. As we expected, the newly 
identified YSOc are well aligned with the "Same" YSOc 
in the diagonal direction, but at fainter ends close to the 
galaxy populated regions. Although in Fig. 4-8, some 
faint YSOc are appeared to be mixed with galaxies, they 
are in fact separated in Multi-D space. Fig. 9 shows 
an example of a 3-D space consisting of IR2, IR4 and 
MP1 magnitudes, in which YSOc and galaxies are better 
separated than in any 2-D space. 

We also compare the infrared luminosity distribution 
of our YSOc list and the "Same" YSOc list (infrared 
luminosity is obtained from J to MP2), and find that 
most of the newly identified YSOc are at the faint end 
of the distribution (Fig. 10a). The distribution of the 
newly identified YSOc is more like an extension of YSOc 
sample at low luminosity end rather than part of the 
galaxy distribution. Therefore, both the YSOc locations 
in MMDs and the luminosity distribution support that 
the Multi-D method are able to identify fainter YSOc 
which are not found by previous works. 

4.1.2. Stellar Surface Density 

Since YSOs tend to form in clusters in molecular clouds 
and background galaxies have relatively random distri- 
bution in the sky, the newly identified YSOc are expected 
to have higher probability to locate in the regions with 
other YSOc if they are real YSOs. We use stellar surface 
density as the parameter to indicate whether a newly 
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Fig. 9. — The 3-D magnitude space with axes of IR2, IR4 and 
MP1. Color points indicate the sources as same as that in Fig. 4 
and all YSOc in five clouds are plotted together. 
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Fig. 10. — (a) The distributions of infrared-luminosity for YSOc 
selected with Multi-D methods (blue), "Same" YSOc (red), YSOc 
with PAH feature (green), and background galaxies (black dashed 
line). The background galaxies are sources with both Gal prob 1 
and Gal prob 2 < 10 and their distribution is scaled by a factor of 
0.1. (b) Same as (a), but only for Class and Class I sources with 
M P2 detections and the black line are VeLLO candidates selected 
by |Dunham et al.|p008} . 



identified YSOc are located within YSO clusters. Stellar 
surface density is defined as 

- 1 



(2) 



where n is the n-th closed star and r n is the dista nce to 
the n-th closed star ( |Gutermuth et al.1 2005, 2009) . Here 
we adopt n=6, same as that used in (Gutermuth et al. 



2009 1 , as a surface density reference. The stellar surface 
density for each source is calculated from where they are 
with respect to the "Same" YSOc population, and the 
results are listed in column 7 of Table 2. In Fig. 11, 
we compare the stellar surface density for the "Same" 
YSOc sample, newly identified YSOc, and background 
galaxies for each cloud. We select background galaxies 
as sources with both Gal prob 1 and Gal prob 2 larger 
than 10. We use Kolmogorov-Smirnov test to determine 
whether the stellar surface density distributions of our 
newly identified YSOc is similar to that of "Same" YSOc 
sample or galaxies. 

The results (P-value) of Kolmogorov-Smirnov test are 
shown in Table 3. If the P-value is significantly larger 
than 0.05, we cannot reject the hypothesis that the dis- 
tributions of the two samples are the same (SciPy Refer- 
ence Guide, 2012). For Perseus, Serpens and Ophiuchus, 
the P-valucs calculated from newly identified YSOc and 
"Same" YSOc are larger than that from newly identi- 
fied YSOc and galaxies. These results suggest that the 
stellar surface density distributions of our newly identi- 
fied YSOc are like to that of "Same" YSOc rather than 
background galaxies in Perseus, Serpens and Ophiuchus. 
In Chamaeleon II and Lupus, the distributions of both 
newly identified and the "Same" YSOc sample show low 
stellar surface density; therefore, the distributions have 
low statistical significance in determining whether newly 
identified YSOc are similar to real YSOs or background 
galaxies. 

4.1.3. YSO Probability 

If YSOs have similar SEDs, they should distribute con- 
tinuously in a rather concentrated region in Multi-D color 
space (Fig. 12). We examine this property by calcu- 
lating a new parameter, YSO probability, in Multi-D 
color space, which is calculated in the way similar to 
galaxy probability in Multi-D magnitude space. Here 
we calculate the YSO probability from the number dis- 
tribution of our identified YSOc in Multi-D color space 
constructed with four axes (IR1-IR2, IR1-IR3, IR1-IR4 
and IR1-MP1). Similar to the galaxy probability calcu- 
lation, the subarrays (2- and 3-D arrays) are also used 
for YSOc with detections only in three and four of the 
all five bands. The cell size is 0.2 magnitude as same 
as that used in galaxy probability calculations, and the 
a of smoothing beam are 2, -y/3/4 x 2 and a/2/4 x 2 
cells for 4-D, 3-D and 2-D arrays, respectively. The YSO 
probability for each source is listed in column 8 of Ta- 
ble 2. YSOc located closer to the center of the YSOc 
populated region in Multi-D color space will have higher 
YSO probability, and only isolated YSOc will have YSO 
probability equal to 1, which is the minimum value for 
YSO probability. Out of 1313 identified YSOc, only 49 
YSOc have YSO probability = 1 and 18 of them arc 
newly identified YSOc. The chances for these "isolated" 
YSOc in Multi-D color space not being real YSO are 
higher; however, it is still possible that small numbers of 
YSOc contain unusual SED features. Since the SEDs of 
YSOs are not exactly identical, lower YSO probability 
suggests that the type of SEDs are less common if they 
are indeed YSOs. Note that the YSO probability can- 
not be used to identify YSOc, because galaxies occupy 
similar color space. 
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Fig. 11, — The histograms of stellar (YSOc) surface density for sources in all clouds. The histograms indicate the "Same" YSOc (red), 
newly identified YSOc without PAH_em feature (blue) and newly identified YSOc with PAH_em YSOc (green), respectively. Black dashed 
line is the background galaxy found in the direction of the cloud and the number is multiplied by a factor of 0.2. 



TABLE 3 

K-S TEST RESULTS FOR STELLAR SURFACE DENSITY - P VALUE 





Perseus 


Serpens 


Ophiuchus 


Lupus 


Chamaeleon 


newly identified YSOc and "Same" YSOc 


0.88 


0.08 


0.21 


0.01 


0.07 


newly identified YSOc and galaxies 


0.00 


0.00 


0.00 


0.06 


0.37 



TABLE 4 

True YSO missed by our selection method 



R.A. 


Dec 


reason description 


degree 


degree 





193.3217719 -77.1196340 DK Cha. This source saturates at IR2, IPv4 and MP1 bands and is not a HREL source. We consider 
only HREL sources (see §2.1). 

246.6100815 -24.4083265 VLA1623. Since it has only been detected at MP1 and longer wavelength, lack of data excludes it 
from our selection process. 

248.0942748 -24.4755023 IRAS16293-2422. IRAS16293-2422A and B are not resolved at MP1. IRAS16293-2422A is not 
detected in IRAC bands, while IRAS16293-2422B is barely detected from IR2 to IR4. These sources 
are not HREL sources. 
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TABLE 5 

YSOC EXCLUDED FROM THE C2D LISTS 



number 


R.A. 


Dec 


c2d classification* Type b 


reason description 




degree 


degree 








Perseus 


1 


51.9117838 


30.2160953 


YSCXred 


B 


Its MP1 flux is band-filled and its location is 
within the PSF of a "Same" YSOc (51.9128437, 
+30.2175414) at 6 arcsec away. 


2 


52.2383644 


31.2386438 


YSOc_star+dust(MPl) 


B 


Its MP1 flux is band-filled and its location is 










within the PSF of a "Same" YSOc (52.2390096, 
31.2377395) at 4 arcsec away. 


3 


52.2578166 


31.2814344 


YSOc_star+dust(IR2) 


B 


MP1 locates in the diffraction spike of an ex- 










tremely bright source and thus possibly has wrong 
flux measurment. 


4 


52.2693214 


31.3682561 


redl 


C 


This is not a HREL source. It is only detected in 
MP1 and the MP1 flux appears to be affected by 
the cloud emission. 


5 


52.2957986 


31.3072254 


YSOcred 


B 


Its IR2, IR3, IR4, and MP1 flux arc band- 
filled, and its location is between two "Same" 
YSOc (52.2944829, 31.3057251) and (52.2969172, 
31.3087350). 


6 


52.3229300 


31.4634218 


YSOcrcd 


B 


It appears like a jet knot 5 arcsec from a "Same" 
YSOc (52.3215278 31.4629087). 


7 


53.1734652 


31.1789357 


YSOc 


B 


Its MP1 flux is band-filled and its location is 
within the PSF of a "Same" YSOc (51.9128437, 
30.2175414) at 2 arcsec away. 


8 


55.9791799 


32.0175333 


YSOc_star+dust(IR2) 


B 


This is a jet knot of HH211 close to the driven 










source. HH211 has only been detected at MP2 
or longer wavelength, thus is not identified by our 
work. 


9 


55.9854925 


32.0146713 


red 


B 


This is a jet knot of HH211 close to the driven 
source. HH211 has only been detected at MP2 
or longer wavelength, thus is not identified by our 
work. 


10 


55.9886063 


32.0131886 


red 


B 


This is a jet knot of HH211 close to the driven 
source. HH211 has only been detected at MP2 
or longer wavelength, thus is not identified by our 
work. 


11 


55.9901535 


32.0124525 


Gale 


B 


This is a jet knot of HH211 close to the driven 
source. HH211 has only been detected at MP2 
or longer wavelength, thus is not identified by our 
work. 


12 


55.9909120 


32.0534315 


red 


C 


This is not a HREL source. Although this source 
has detections with S/N=3, 2 and > 7 at IR2, IR3 
and MP1, we think it is not a real source due to 
the contamination from cloud emission. 


13 


55.9975507 


32.0099032 


YSOcrcd 


B 


This is a jet knot of HH211 close to the driven 
source. HH211 has only been detected at MP2 
or longer wavelength, thus is not identified by our 
work. 


14 


56.0797337 


32.2883878 


YSOc 


A 


Giant star in this work (§3.3). 


15 


56.1498481 


32.1567564 


YSOcstar+dust(IR2) 


B 


IR4 and MP1 detections seem to be part of the 










clouds. The source appears to be star-like at 
shorter wavelengths. 


16 


56.4473566 


31.7198486 


YSOcstar+dust(IR4) 


C 


This is not a HREL source. 


Lupus 


17 


234.7014819 


-34.6772876 


YSOc_PAH-em 


B 


This source is elongated in IR1 image, thus we 
identified it as a galaxy. 


18 


242.2427972 


-39.1265213 


YSOc 


A 


Giant star in this work (§3.3). 


19 


242.3921400 


-39.2283528 


YSOc 


A 


Giant star in this work (§3.3). 


Ophiuchus 


20 


246.5609716 


-24.4187556 


redl 


U 


This is not a HREL source. It is only detected 
in MP1 and we think the MP1 flux is contributed 
from cloud. 


21 


246.7462569 


-24.5842515 


red 


B 


Its MP1 flux is band-filled and its location is within 
the PSF of a newly identified YSOc (246.7464700, 
-24.5829609) at 6 arcsec away. 


22 


246.8052355 


-24.6926085 


YSOc 


A 


Giant star in this work (Gaint star GY 232 iden- 
tified by Luhman et al. 1999). 


23 


246.8747631 


-24.5601131 


YSOc 


A 


Giant star in this work (§3.3). 


24 


246.9323892 


-24.7188227 


YSOc 


A 


Giant star in this work (§3.3). 
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TABLE 5 

YSOC EXCLUDED FROM THE C2D LISTS 



number R.A. Dec c2d classification 51 Type b reason description 

degree degree 



Serpens 



9^ 


977 


U40OUUU 


-U.UZ / OUU'i 


VCJO/-' e+«r-L-rliic-r71\/TP1 "l 
I D WC_SId,r-)-QUht ^lVl_r 1 J 


A 


Giant star identified by Olivcira ct al. (2009) 


26 


277 


1141023 


-0.1972144 


YSOc_star+dust(MPl) 


A 


Giant star identified by Oliveira et al. (2009) 


27 


277 


2270870 


0.4812201 


YSOc 


A 


Giant star in this work (§3.3). 


28 


277 


2837039 


-0.1269813 


YSOc_star+dust(MPl) 


A 


Giant star identified by Oliveira et al. (2009) 


29 


277 


2876881 


0.5244498 


Gale 


B 


This is a jet knot of "Same" source (277.2877868, 
0.5256632) close to the driven source. 


30 


277 


3855879 


-0.2231446 


YSOc_star+dust(MPl) 


A 


Giant star identified by Oliveira et al. (2009) 


31 


277 


5219079 


0.6846068 


red 


E 


This source locates in the edge of the mapping area 
of c2d, thus not observed in IR1 and IR3. We do 
not consider such sources in our analysis. 


.32 


277 


5237366 


0.6588334 


rising 


E 


This source locates in the edge of the mapping area 
of c2d, thus not observed in IR1 and IR3. We do 
not consider such sources in our analysis. 


33 


277 


5451572 


0.7835665 


YSOc_star+dust(IR4) 


B 


We identify this source as a galaxy due to its 
galaxy probability (Gal probl = 1.07 and Gal 
prob2 =1.04). 



a If the c2d source type does not contains "YSO", the source is added to the YSOc list by |Evans et al.| | |2009t based on ancillary data 
from literature. 



Reasons for excluding the source from our YSOc list. 
A: This source is classified as a giant star by either this work or Oliveria et al. (2009). 
B: This source is excluded in image checking process. 

C: This source is not considered in our identification process. It has detections at fewer than three bands or is not a HREL source 
(see §2.1). 

D: This source has both galaxy probabilities, Gal prob 1 and Gal prob 2, larger than 1. 
E: This source is not observed at IR1 and IR3 because it is in the edge of c2d map 
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Fig. 12. — Color-color diagram using IR1, IR2, IR3 and IR4 
bands. Color points indicate the sources as same as that in Fig. 4 
and all YSOc in five clouds are plotted together. 



4.2. Missed or Excluded YSOc from The c2d Lists 

In this section, we discuss the c2d identified YSOc 
but missed through our selection processes. There are 
36 sources missed in total. We recover only three of 
them back to our YSOc candidate list (Table 2) because 
they are previously identified well-known sources. These 
three sources are not identified in our identification pro- 
cess because they are not HREL sources (Table 4); DK 
Cha is saturated at several Spitzer bands, and VLA1623 
and IRAS 16293 are only detected at MP1. For the re- 
maining 33 sources, they ar e separated into two g roups: 
23 sources are identified by |Harvey et al. ( |2007a[ ) using 



several CMD criteria ("c2d classification" starting with 
"YSOc" i n column 4 of Tab le 5), and 10 sources are 
added by Evans et al. (2009) with additional observa- 



tions from othe r literat ures ( iDunham et al. || 2008 l |Enoch 
et al.||2009| |j0rgensen et al.)|2007| |2UU8|) . Here we justify 
why they are not likely to be true YSOs, and the exact 
reasons for the removal of each source are listed in Table 
5. 

First , for the 23 YSOc identified by |Harvey et al. 



(2007a 



1. 



Eleven sources may be giant stars: seven sources are 
classified as giant stars by this work, in cluding a giant 
star G Y232 independently identified by |Luhman et al.] 
(1999 ) with spectrosc opic data, and four s ources are 
identified as giants by Oliveira et al. (2009). 



2. Ten sources are removed by image checking process: 
seven sources have flux contamination from nearby 
sources or cloud emission, one source is elongated like 
a galaxy, and two sources are jet knots from HH211. 

3. One source is not a HREL object. It is also very dis- 
tant from the cluster regions with low stellar surface 
density in Perseus and slightly extend, thus we suggest 
it is a background galaxy. 

4. One source has Gal Prob >1 (1.07). 



Seco nd, for the ten c2d YSOc added by |Evans et al. 
( 2009 1 but removed by our work, 



Three sources were identified as jet knots in HH211 
system by image checking. Note that HH211 itself 
also cannot be identified through our procedure, since 
it has only been detected at MP2 and longer wave- 
lengths. 
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2. Two sources locate very close to selected YSOc, and 
their flux appear to be contaminated by the nearby 
YSOc. 

3. Five sources have not been considered in our identifi- 
cation process. Three of them are not HREL sources, 
and they all appear to be extended and seems to be 
part of the cloud structure. The remaining two are 
located at the edge of survey regions without IR1 and 
IR3 observations, which are not considered in our anal- 
ysis. 

4.3. PAH Emission 

From SWIRE data, we find that the galaxies with a 
peak at IR4 in SEDs (PAH emission) are usually bright, 
thus are much difficult to be eliminated. The PAH emis- 



sio n are usually seen in star- forming galaxies ( Evans et 



al. |2009 ), thus the c2d project label such sources with 
TAH_em" (col. 9 in Table 2) usin g the color criteri a 
[3.6]-[4.5] > 1.5 and [5.8]-[8.0] < 0.6 ( |Evans et al.]2007| . 
There are 37 sources with "PAH_em" label out of our 
1313 YSOc (30 from newly identified YSOc and 7 are 
"Same" YSOc). The YSOc satisfying the PAH-em con- 
ditions are not excluded in neither the c2d YSOc list 
nor our YSOc list, because they are not clearly identified 
as galaxies. More observations toward these sources are 
necessary for identifying their nature, and we should use 
these sources with caution. Although those YSOc could 
be wrongly identified, the small fraction of 2.8% should 
not significantly influence our statistical studies. 

4.4. Uncertainty From The Selections of Bands for 
The Two 5-D arrays 

Because we use two specific 5-D arrays to identify 
YSOc instead of a 10-D array, it can result in an un- 
certainty in our YSOc sample from the specific selec- 
tion of bands. To understand the effect, we construct 
another three Multi-D arrays to select YSOc and com- 
pare the results with our YSOc identified using bands 
of J, K, IR2, IR4 and MP1 (1277 sources). These three 
Multi-D arrays are constructed using bands of A=(J, H, 
IR1, IR3 and MP1), B=(J, Ks, IR1, IR4 and MP1) and 
C=(J, IR1, IR2, IR4 and MP1), respectively. They all in- 
clude J and MP1 bands to represent the SED structures 
with the longest coverage of wavelengths and subarrays 
are also constructed for each main array. Following the 
whole selection process, there are 7, 7 and 2 additional 
sources selected using A, B and C arrays, respectively. 
One source is selected by all three spaces and another 
one is selected by both A and B spaces; thus, totally 13 
new sources are selected using these three spaces. These 
13 sources all have galaxy probability greater than 0.4 
which is corresponded to ~ 2 cells (0.4 magnitude) from 
a boundary of galaxy populated region with our Gaus- 
sian beam. Therefore, we suggest that the specific band 
selection could result in < 1% difference compared to our 
original YSOc selection and the difference only occurs in 
the region near the boundaries between YSOc and galax- 
ies. 

5. DISCUSSION 
5.1. New YSO Identification Method and Results 



We develop Multi-D method for separating two kinds 
of sources which have similar SEDs and are indistin- 
guishable in CMDs and CCDs. In this paper, we use 
Multi-D method to identify YSOc in molecular clouds 
with Spitzer data by comparing to galaxy sample from 
SWIRE. A large number of new YSOc are identified 
and our analysis indicates that they are much like faint 
YSOs not selected before (Fig. 10). This result suggests 
that Multi-D method is very powerful in identifying faint 
YSOs that have been detected. However, there are sev- 
eral caveats in Multi-D method. 

1. The SWIRE data may not represent a complete col- 
lection of galaxy SEDs. In addition, in the direction 
of the observed clouds, the brightest galaxies may be 
brighter than all SWIRE galaxies due to their prox- 
imity. Thus, if SWIRE data are incomplete galaxy 
samples, some true galaxies will not be removed and 
will be identified as YSOc. 

2. Some YSOs may have similar SEDs and brightness to 
those of galaxies. If a YSO has SED and brightness 
similar to galaxies in the comparison galaxy sample in 
the analysis wavelengths, it will have a galaxy proba- 
bility larger than one and will be removed from YSOc 
list. We found that our YSOc and galaxy sample are 
continuously distributed in MMDs and 3-D plotting 
(Fig. 4-9), which implies that there may be faint 
YSOs located in the galaxy populated region and can 
not be identified with Spitzer data only. However, no 
available methods can identify such sources with only 
photometry data. 

3. Some sources may have high uncertainty in photom- 
etry which will affect the YSO identification process. 
Sources with wrong flux measurements will have pe- 
culiar SED structures, which makes them locate far 
from the galaxy populated region in Multi-D space 
and thus be classified as YSOc. We have designed an 
image checking process to reduce such cases (§3.4). 

Although Multi-D method has these caveats, such prob- 
lems are unavoidable and also exist in those methods us- 
ing CCD and CMD criteria. Multi-D method does allow 
us to obtain the most complete YSOc list from provided 
data set. 

5.2. Newly Identified Very Low Luminosity Object 
(VeLLO) Candidates 

VeLLOs are faint embedded protostars (L int < 0.1 L ), 
therefore they are difficult to be identified. Recent works 
suggest that VeLLOs a re likely to be protos tars at qui- 
escent accreting phase (Dunham et al. 20101. However, 
this does not exclude the possibility that some VeLLOs 
could be faint protostars at early evolutionary stages or 
even "first cores" which is the transition phase between 
starless cores and Class sources. Since our method dis- 
covers more faint protostars compared to previous works, 
we here examine whether we can identify more VeLLOs. 
The most thorough su rvey of VeLLOs is carried out 



by Dunham et al. (20081 using c2d data. Several other 
sources are claimed to be VeLLOs in the re cent years 

Since a 



(Pineda et al. 2011 



Kauffmann et ah] 2011). 
large number of taint YSOc are identified using Multi- 
D method in this work, we are in good position to add 
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TABLE 6 
VeLLO candidates 



R.A. 


Dec 




t 70/im 
int 

L Q 


c2d classification 


Ay 


degree 


degree 


Lq 




mag 



52.2752091 
52.3348374 
52.7512777 
53.2410000 
56.3075775 
194.2566223 
277.2287724 



30.5108863 
31.2139919 
30.9369012 
31.1022931 
32.2027786 
-76.8097467 
0.3090777 



0.017 (0.007) 
0.016 (0.006) 
0.009 (0.004) 
0.010 (0.004) 
0.050 (0.02) 
0.028 (0.007) 
0.034 (0.003) 



0.023 (0.006) 
0.034 (0.008) 
0.011 (0.003) 
0.016 (0.004) 
0.054 (0.013) 
0.049 (0.008) 
0.038 (0.005) 



YSOc_star+dust(IR2) 



YSOcred 
YSOc_star+dust(IR2) 

YSOc_red 



2.4 
9.5 
0.9 
13.7 
3.0 
2.4 
15.9 



1 




! 3 


4 


i 5 


i 6 




10° 10 1 10 2 10° 10 1 10 2 
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Fig. 13. — The SEDs of VeLLO candidates identified from our 
YSOc list (see §5.2 for selection criteria). Number 3 is newly identi- 
fied from our work, and the rest of VeLLOs have been identified by 
Dunham et al. (2008) but was classified as less confident VeLLOs. 



reliable VeLLO candidates to the collection. Fig. 10b 
shows that some newly identified Class and Class I 
YSOc (with MP2 detections) have lower infrared lumi- 
nosity (luminosity between 1.25-70/im) than that of the 
mo re reliable VeLLO can d idates (group 1-3) ide ntified 



p 1 



Dunham et al. (2008) use a 



by Dunham et al 
relation to translate the M 
nosity Lj nt = 3.3 x 10 
in cgs units (ergs cm~ 2 s _1 ), based on the SED model 
from Monte Carlo dust radia tive transfer code RADMC 
( Dullemond & Dominik|2004 ) . We adopt this relation to 



(|2008|). 

2 flux to the internal lumi- 
8 F^ 94 L where F 70 is MP2 flux 



estimate L; nt tor our Class and I candidates and find 
32 sourc es with Lj nt < .1 L^ in addition to VeLLOs 
found in Dunham et al. (20081. In order to make sure 



the s elected VeLLOs are yo ung YSOs, two more criteria 
from |Dunham et al. ( 2008 ) are adopted to ensure that 
the SEDs rise (1) from the longest detected IRAC wave- 
length to MP1 and (2) from MP1 to MP2, which reduce 
the number of VeLLO candidates to seven (Table 6). 

Their SEDs are shown in Fig. 13. Among these seven 
candidates, only one is newly identified by us. The 



oth er six cadidates have been identified by Dunham et| 
al. (2008), but they are in less confident group 4-7 corre- 
sponding to sources not obviously associated with high 



column density regions. Since one possible nature of VeL- 
LOs is proto brown dwarfs which could be ejected from 
their parent cores during formation processes, VeLLOs 
may not necessarily be embedded in high density regions. 
Therefore, we do not exclude these sources from our 
VeLLO candidate list. In short, with Multi-D method, 
faint protostar candidates can be selected naturally. 

5.3. Statistical Properties 

The significant increase of YSOc numbers in molecular 
clouds can alter statistical properties, such as SFR and 
lifetimes of YSOs in different evolutionary stages. The 
variations of these two properties are analyzed in this 
section. 

5.3.1. Comparing Star Forming Rate with Theoretical 

Models 

The accuracy of SFR is critical in distinguishing 
whether the star formation in GMCs is dominated by tur- 
bulence or magnetic fields. Without any supports, typi- 
cal GMCs should collapse on its free-fall tim e scale, which 
results in a S FR of roughly 250 M yr _1 (Krumholz & 
McKee||2005| ). Howeve r the SFR in Milky Way is mea- 
sured as ~ 3 M yr" 1 ( |McKee fc Williams||l997| ) . Thus, 
supporting forces, such as turbulence and/or magnetic 
fields, are required to reconcile the difference between 
theory and obse rvations. To study how the turbulence 
affects the SFR iKrumholz & McKee] ( [2005[) defines the 
dimensionless SFR per unit tree-tall time, SFR ff , 



SFR ff = 



SFR x i ff 
M (cloud) 



(3) 



where tff is the free fall time of the cloud and M(cloud) is 
the mass of the cloud, and derive an analytic expression 
for SFRff in a supersonic turbulent median, 



SFRff 



0.014(^)-°- 68 (^)-°- 32 , 



(4) 



where a 



is the virial parameter defined by |Bertoldi fc| 



McKee ( 1992 ) measuring the ratio of kinetic energy and 
gravitational energy of a clumpy, and A4 is the Mach 
number. The definition of a V i r is a V i r —5a^ ot R/ (GM) 
where o~tot is velocity dispersion from thermal and tur- 
bulent velocities over entire cloud, R and M are radius 
and mass of the cloud, respectively. The turbulence dom- 
inated model suggests that a small a V i r can result in a 
large SFRff. However, magnetic field domina ted models 
predict that the SFRff is proportional to a v 
fc Tan|[2007l ), which is 



( Krumholz 



SFRff 0.01a t 



(5) 



1G 



TABLE 7 
Facts and parameter, about clouds 



Cloud 




N(YSO) b 


M (cloud) a 




Av a 


SFRff 






pc 


# 


M© 


pc 2 


km s -1 






Perseus 


250 (50) 


452 


4814 (1925) 


62.0 (24.8) 


1.54 (0.11) 


0.050 (0.034) 


0.46 (0.21) 


Serpens 


260 (50) d 


295 


2016 (775) 


17.5 (6.7) c 


2.16 (0.01) 


0.047 (0.030) 


1.14 (0.49) 


Ophiuchus 


125 (25) 


343 


2182 (873) 


27.3 (10.9) 


0.94 (0.11) 


0.068 (0.045) 


0.25 (0.13) 


Lupus I 


150 (20) 


20 


250 (58) 


4.5 (1.2) 


1.9 (0.19) f 


0.026 (0.012) 


3.62 (1.30) 


Lupus III 


200 (20) 


60 


443 (102) 


8.5 (1.7) 


1.7 (0.17) f 


0.054 (0.018) 


2.25 (0.68) 


Lupus IV 


150 (20) 


13 


119 (28) 


1.7 (0.5) 


1.7 (0.17) f ' 


0.025 (0.011) 


3.74 (1.34) 


Chamaeleon II 


178 (18) 


27 


426 (86) 


7.4 (1.5) 


1.2 (0.4)s 


0.023 (0.008) 


1.09 (0.77) 


a The area and Av 


of clouds are 


from|Evans et al.|l|2009[. 











c The area is calculated from extinction map with Av > 2 magnitude area. 

_ We use th e dist ance uncertain t y of S e rpens with a wide range to cover the distance measurement mention inlEvansI 
|et al.|p009} , i.e., |Straizys et al.| jl996[ ; |Eiroa et al.| l |2008[ 

The extinction map of Serpens has all pixel with Ay > 2 magnitude, thus the area is from Spitzer observing region. 



^ The Av in Lupus is from pTachihara et al 



i 19961 



s The Av in Chamaeleon is from [Vuas-^oaset"al.| j!994[ 



The parameter values Krumholz & Tan (20071 used to 
derive this relation could vary more than an order of 
magnitude, resulting in a wide range of possible value for 
the coefficient. Nevertheless, the linear proportionality of 
SFRff and a V i r should be unaffected. Hence, comparing 
the relation of SFRff and a V i r may provide us the hint 
that which mechanism dominates in molecular clouds. 

The c2d data and our new YSOc catalogs provide nec- 
essary information to compare SFRff and a V i r . S FRff is 
estimated f ollowing the same analysis method in Evans 
et al. ( 2009 1 . Assuming that the mean mass of a YSO M* 

2Myr, 

(G) 

(7) 



U75 Mq and the period of star formation T 
the star formation rate can be estimated from 



SFR 



N(YSO) x M* 
T 



and 



iff = 34Myr/Vh, 



where N(YSO), M*, T and n are the number of YSO, 
mean mass of YSO, period of star formation and number 
density of cloud, respectively. Inserting Equation (6) and 
(7) into Equation (3), we are able to obtain the SFRff 
from N(YSO), M(cloud) and n. We calculated n from 
the cloud mass and surface area by assuming a spherical 
cloud and the cloud mass was obtained from Table 1 in 
The cloud mass and the surface area 



Evans et al](|2009) 

UC1! 



of the clouds are both obtained from the c2d extinction 
maps in area with Ay > 2. We use the mass of the clouds 



from the table 1 in Evans et al 
surface area of the clouds 



( 2009 1 and calculated the 
which are both from the c2d 
> 2 region. The errors 



extinction maps in area with Ay 
of SFRff are from error propagating and only uncertainty 
considered here is the uncertainty in cloud distances. 

We calculate a v i r using a V i r =5a^ ot R/ (GM) . R and M 
are the same to that used in SFRff calculations. The 
Utot are velocity dispersions from 13 CO J = 1 — > 
observations and are averaged over a full map of the 
cloud. Fo r Perseus, Serpens a nd Ophiucuhus, we use the 
a tot from"" 
from the ' 



Evans et al. 



( 2009 ) which obtained the values 
he Coordinated Molecular Probe Line Extinc- 
tion Thermal Emission Survey of Star Fo rming Regions 
(COMPLETE) project ( |Ridge et al. 1 [2006|). For Lupus 
and Chamaeleon, we use the a to t from |Tactnhara et al.| 




o.i 



0.05 



0.02 




Lupus I 



LSipis IV - 



Fig. 14. — The plot of a v i r to SFRff. The solid line and dashed 
line represents the best fitting curve of eq. 4 (turbulence model) 
and eq. 5 (magnetic field model), respectively. 



(19961 and Vilas-Boas et al. (19941, respectively, which 



also use 13 CO J = 1 — > as the tracer. The calculated 
a v i r for all clouds are shown in Table. 7 and the propa- 
gating errors are from the errors of distance and a to t in 
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the cloud. 

Fig. 14 shows the relations between a V i r and SFRff 
for all clouds in this paper. We fit the data with the 
predictions for the turbulence model and the magnetic 
field mo del, i.e., Equation (4) a nd (5), respectively, de- 
rived by Krumholz & Tan (2007). For the magnetic field 
model, we set the coefficient to be a variable in Equation 
(5). The reduced x 2 are 3.0 and 5.6 for the turbulence 
model and the magnetic field model, respectively, which 
implies that the data are more consistent with the turbu- 
lence model than the magnetic field model. Our result 
hints that turbulence dominates the star formation in 
large scales such as GMCs. In addition, in Fig. 14, the 
SFRff of Lupus are all larger than that of the best fitting 
curve of the turbulence model, especially for Lupus III. 
This could result from our assumption of a mean mass 
0.5 M Q for YSOc in all clouds, but t he mean stellar mas s 



Class I/O 



in Lupus may be closer to 0.2 M Q ( |Merm et al. 2008). 
The mean stellar mass is 0.6 or 0.8 Mq for different evo- 



lutionary tracks in Serpens ( Oliveira et a l. 2009) and i s 
0.52 ± 0.11 M Q in Chamaeleon 11 ( |Spezzi et al"||2008 l. 
Better estimates of the mean stellar mass in all clouds 
will provide a more accurate SFRff and give a stronger 
conclusion. 




Fig. 15. — The populations of all sources with age indicator, a, 
in the five clouds. Blue and red lines indicate our YSOc and the 
"Same" YSOc, respectively. 



5.3.2. Lifetime in Different Evolutionary Stages 

The lifetimes of YSOs in different evolutionary stages 
are estimated from the numbers of YSOs in each stages, 
and thus the variation of YSOc sample may result in dif- 
ferent lifetime estimations. YSOs are commonly classi- 
fied into four evolutionary stages Class 0/1, Flat, Class II 
and Class III from young to old using the spectral slope, 
a, which is the best fit slope from K to MP1 band of the 
SED. The definition of a and the classification criteria 



are from Greene et al. ( 1994 ) 



d log(AS(A)) 
d log(A) 



(8) 



and 

Class 0/1: 0.3 < a 
Flat: -0.3 < a < 0.3 
Class II: -1.6 < a < - 
Class III: a < -1.6 



0.3 



Although it has been demonstrated that the SED mor 
phologies of YSOs are affected by geometry such as in- 



al. 2003a 



still 



Robitaille et al. 


2006 


use a as a roug 


h age 



et 



_ ., fiapsi et al.||2008 | ,, 
i age indicator as there is 
no better alternatives. Our large YSO sample size may 
mitigate the geometry effect. 

Fig. 15 shows the distributions of a in each clouds 
and the fractions of YSOc in different evolutionary 
stages, and th e latter numbers are almost the same to 
the results in Evans et al. (2009) in either individual 
cloud or the whole data set, implying that the estimates 
of lifetimes of each evolutionary stages are not changed 
(Table 8). Fig. 15 also shows that the distributions 
of a in different clouds are different, which indicates 
that those clouds may be in different evolutionary 
stages. Because Perseus molecular cloud consist of two 
major YSO clusters NGC1333 and IC348, we separated 
it into two regions, western and eastern Perseus, by 



R.A.=54.3°. Our results suggest that western Perseus 
(NGC1333) has the highest fraction of Class 0/1 YSOc, 
while Lupus has the highest fraction of Class III YSOc. 
Therefore, western Perseus is the youngest cloud and 
Lupus is the oldest cloud among c2d-surveyed regions, 
while the other clouds are at similar evolutionary stages. 



6. SUMMARY 

We have developed "Multi-D method" to identify 
YSOc from star-forming regions with reliable photom- 
etry measurements from multiple bands. Main-sequence 
and Giant stars are first eliminated from the observed 
sources, and the rest of the sources are compared to a 
galaxy sample, such as Spizers SWIRE dataset, in multi- 
dimensional magnitude space. We demonstrate that this 
Multi-D method which uses multi-band photometry si- 
multaneously can identify all possible YSOc above the 
galaxy confusion limit and recover those could be missed 
by using certain set of color-magnitude or color-color cri- 
teria. We identify 1313 YSOc from Spitzer's c2d high re- 
liability (HREL) catalogs, which is 28% m ore than what 
have been identified by the c2d project (Evans et al. 
20091. The increase amount of YSOc suggests the fol- 
lowing results: 



1. The increase amount of YSOc directly increases the 
SFR estimated in each region by the same percent- 
age. We further compare the relation between a V i r 
and SFRff to theoretical models, and our results are 
more consistent with the prediction of supersonic tur- 
bulence dominated model than that of magnetic field 
dominated model. 

2. Although we identify 28% more YSOc than those 
listed in the c2d catalogs, the lifetimes of YSOs in 
different evolutionary stages are unchanged since the 
fractions of YSOc in different Classes are almost the 
same. 
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TABLE 8 

Numbers of YSOc by clouds and class 



Cloud 




This work 






jEvans et al. 2009 1 




1/0 


Flat 


11 


III 


1/0 


Mat 


11 


III 


Perseus 


99 


49 


272 


49 


87 


42 


225 


31 


Serpens 


55 


36 


170 


35 


36 


23 


140 


28 


Ophiuchus 


72 


60 


191 


58 


35 


47 


176 


34 


Lupus 


12 


12 


64 


41 


5 


10 


52 


27 


Chamaeleon II 


7 


3 


21 


7 


2 


1 


19 


4 


Total 


245 (19%) 


160 (12%) 


718 (55%) 


190 (14%) 


165 (16%) 


123 (12%) 


612 (60%) 


124 (12%) 



3. Our Multi-D method allows us to reliably identify 
more faint YSOc, which can be used to select VeL- 
LOs. Using the 70/im flux to the internal luminosity 



relation suggested by Dunham et al. (20081, we are 
able to find 7 new VeLLU candidates. 
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